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Abstract. We explore and compare the performances of two nonlinear correction 
and scale-dependent biasing models for the extraction of cosmological information from 
galaxy power spectrum data, especially in the context of beyond-ACDM cosmologies. 
The first model is the well known Q model, first applied in the analysis of 2dFGRS 
data. The second, the P model, is inspired by the halo model, in which nonlinear 
evolution and scale-dependent biasing are encapsulated in a single non-Poisson shot 
noise term. We find that while both models perform equally well in providing adequate 
correction for a range of galaxy clustering data in standard ACDM cosmology and 
in extensions with massive neutrinos, the Q model can give unphysical results in 
cosmologies containing a subdominant free-streaming dark matter whose temperature 
depends on the particle mass, e.g., relic thermal axions, unless a suitable prior is 
imposed on the correction parameter. This last case also exposes the danger of 
analytic marginalisation, a technique sometimes used in the marginalisation of nuisance 
parameters. In contrast, the P model suffers no undesirable effects, and is the 
recommended nonlinear correction model also because of its physical transparency. 
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1. Introduction 

The past decade saw an explosion of precision cosmological measurements. Foremost 
amongst these is the observation of temperature and polarisation fluctuations in the 
cosmic microwave background (CMB) by a range of experiments [1-6]. The distribution 
of large-scale structures (LSS) has also been mapped to unprecedented breadths and 
depths by galaxy redshift surveys such as the Two-Degree Field Galaxy Redshift 
Survey (2dFGRS) [7] and the Sloan Digital Sky Survey (SDSS) [8,9]. Together with 
observations of distant type la supernovae [10, 11], these measurements have fostered the 
emergence of a benchmark framework — the adiabatic, nearly scale-invariant, "vanilla" 
ACDM model — based on which one can test for evidence of new physics. 

Clustering statistics of galaxies as probed by surveys like 2dF and SDSS are 
particularly well suited for the exploration of physics that introduce new effects on 
length scales (9(10) — ► (9(100) h^ 1 Mpc. A classic example is the possibility to detect a 
subdominant component of free-streaming hot dark matter (HDM) [12], notably massive 
neutrinos [13,14] and variants such as thermal axions [15-19], light gravitinos [20], 
et cetera. The sensitivity of these surveys at small length scales also lends a greater lever 
arm to the search for features in the primordial density perturbation power spectrum, 
possible remnants of inflationary physics [21,22]. Last but not least, since the power 
spectrum of large-scale structures probes uniquely the parameter combination fl m h, it 
helps to lift the degeneracy between — and hence tighten the constraints on — the physical 
matter density Vt m h 2 and the Hubble parameter h from CMB observations alone. 

Usage of data from galaxy clustering surveys is based on the premise that one 
can reliably predict the distribution of galaxies, at least on a statistical basis, from 
theory. This is complicated by a number of factors. First, galaxies are necessarily 
collapsed objects, i.e., they have undergone a phase of nonlinear evolution. Using them 
as tracers of the underlying matter field implicitly assumes we know how to relate 
the two distributions to one another. A reasonable assumption is that on sufficiently 
large scales the power spectra of galaxies and of the matter field are identical up to 
a constant normalisation, or bias, factor. But this bias relation is expected to become 
scale-dependent when the dimesionless power spectrum of the galaxies A ga i = /c 3 P gal /27r 2 
exceeds unity [23,24]. Indeed, the apparent tension between the 2dF and the SDSS 
galaxy power spectra is now believed to have originated from a more strongly scale- 
dependent bias factor for the red galaxies dominating the SDSS galaxy catalogue [25, 
26]. On the theoretical front, a good deal of recent effort has also been devoted to 
understanding the origin of scale-dependent biasing (e.g., [27,28]). 

Second, galaxy positions are inferred from their redshifts. However, the peculiar 
motions of the galaxies, particularly when amplified by virialisation, induce additional 
Doppler shifts that can potentially obscure the inference. This is known as redshift space 
distortion, and on small length scales requires corrections beyond linear perturbation 
theory. Third, the clustering of the underlying dark matter field itself becomes nonlinear 
on scales k > 0.2 h Mpc -1 . Thus, how reliably one can extract cosmological information 
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from galaxy clustering statistics depends crucially on how well one can model these three 
nonlinear effects. 

While all three effects can in principle be modelled by numerical simulations, these 
simulations, and indeed our understanding of galaxy formation, are not yet at a stage 
where one can reliably predict the power spectrum of galaxies as a function of galaxy 
type and redshift given some underlying cosmological model. In the meantime, nonlinear 
evolution and scale-dependent biasing must be modelled empirically as a systematic 
effect and the associated nuisance parameters marginalised when extracting cosmological 
information from galaxy clustering surveys, especially in beyond-ACDM cosmologies. 

In this connection, Cole et al. [7] recently proposed a correction formula which 
maps directly between the matter power spectrum calculated from linear perturbation 
theory Pu n (k) and the power spectrum expected for the galaxies P ga x(k), 

P ^ {k) = h2 T^ PUk) - (L1) 

The formula is partially calibrated against ACDM-based semi-analytic galaxy formation 
simulations (A n \ = 1.4 for redshift space, and 1.7 for real space), and contains two free 
parameters (b and Q n \) to be fixed by observational data. Equation (1.1) has been 
applied to the galaxy power spectra of 2dF [7] and the SDSS luminous red galaxy 
(LRG) sample [8] to test standard vanilla cosmology. However, there is a priori no 
guarantee that its usefulness extends also to cosmologies beyond ACDM. Indeed, as 
we shall show below, equation (1.1) can be highly pathological when applied to certain 
classes of cosmological models containing a subdominant component of free-streaming 
dark matter. 

In comparison, a conceptually more appealing framework in which to discuss 
nonlinear corrections is the halo model [29-32]. Building on the assumptions of 
(i) hierarchical clustering, and (ii) that galaxies form only inside dark matter halos, halo 
model-based nonlinear corrections can in principle be made applicable to all hierarchical 
CDM cosmologies. The minimal model proposed in references [33-35], for example, 

Pgal(fc) = b 2 P hn (k) + P shot , (1.2) 

where b and P s hot are free parameters, does not demand vanilla ACDM as the sole input 
cosmology. 

In the present work, we explore and compare the performance of these nonlinear 
correction models in some detail. We confront them with various observed galaxy power 
spectra, especially in the context of beyond-ACDM cosmologies. We discuss parameter 
degeneracies and the role of priors on the nuisance parameters. 

The paper is organised as follows. We describe first in section 2 the galaxy clustering 
data sets used in the analysis. Section 3 contains a more detailed discussion of the two 
nonlinear models we wish to explore. In sections 4, 5 and 6 we test the nonlinear models 
against data for standard vanilla cosmology, vanilla with massive neutrinos, and vanilla 
with thermal axions respectively. We conclude in section 7. 
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2. Data sets 

We use the publicly available galaxy power spectra from the following galaxy catalogues. 

2dF This data set comes from the final data release of the Two-Degree Field Galaxy 
Redshift Survey [7]. We use up to 36 data bands, corresponding to redshift space power 
spectrum data for wavenumbers 0.02 < k/h Mpc -1 < 0.18. 

SDSS-2 main Real space power spectrum of the main galaxy sample from the 
Sloan Digital Sky Survey data release 2 [36]. We use up to 19 data bands, i.e., 
0.016 < k/h Mpc" 1 < 0.2. 

SDSS-4 LRG Real space power spectrum of the luminous red galaxies from the Sloan 
Digital Sky Survey data release 4 [8]. The 20 data bands correspond to wavenumbers 
0.012 < k/h Mpc- 1 < 0.2. 

WMAP-3 We also use at times CMB data from the Wilkinson Microwave Anisotropy 
Probe experiment after three years of observations [37-39], mainly for the construction of 
priors on certain cosmological parameters. This calculation is performed using version 2 
of the likelihood package provided by the WMAP team on the LAMBDA web page [40]. 

3. Two nonlinear models 

3.1. The Q model 

We refer to the correction formula (1.1) as the Q model. For ACDM cosmologies, galaxy 
formation simulations suggest that the parameter A n \ can be held fixed at A n \ = 1.4 in 
redshift space and at A n \ = 1.7 in real space [7]. The parameter Q n \, however, exhibits 
a strong dependence on the galaxy type. Fitting the Q model to the 2dF galaxy power 
spectrum at 0.02 < k/h Mpc -1 < 0.3, Cole et al. found Q n \ = 4.6 ±1.5 for vanilla 
cosmology [7]. Tegmark et al. [8] applied the same model to SDSS-4 LRG, for which 
Q n \ = 30 ± 4 provides a good fit 4 

Note that the case of Q Q \ = is not equivalent to no nonlinear correction, since 
a nonzero A Q \ parameter also modulates the power spectrum in a scale-dependent way. 
For Q n \ < 7, the Q model suppresses the power spectrum at k < 0.2 h Mpc -1 . However, 
for Q n \ values as large as 30 such as required by the SDSS-4 LRG, the main role of the 
Q model is to add power at k > 0.07 h Mpc -1 . 

| Tegmark et al. [8] used A n \ = 1.4, although strictly speaking A n \ = 1.7 is the more appropriate 
value for the real space power spectrum of SDSS-4 LRG. However, at the present level of precision, our 
tests show that adopting the correct ^4 n i = 1.7 only causes a statistically insignificant upward shift in 
the best- fit Q n \ (AQ n i ~ 2), while the estimates of other cosmological parameters remain unaffected. 
Henceforth we shall use exclusively A n \ = 1.4 for both real and redshift space power spectra. 
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3.2. The P model 

The functional form of the Q model (1.1) was recently criticised in reference [28] for 
its incorrect dependence on k. Specifically, it lacks a constant term to account for 
the presence of non-Poisson shot noise, a generic consequence of the assumption that 
galaxies form exclusively in halos. Indeed, from the halo model [29-32], one should 
expect a correction formula with the skeletal form, 

P ga ,i{k) = P 2h {k) + P lh (k). (3.1) 

Here, the two-halo term, P2h{k) = b 2 (k)P\i n (k), arises from correlations between galaxies 
in two different halos, and approximates the familiar linear bias relation on large scales, 

P2M » b 2 P hn (k). (3.2) 

The one-halo term 

Pih(fc) ~ const. = P shot (3.3) 

accounts for correlations within the same halo and is approximately independent of 
the exact spatial distribution of the galaxies within a halo provided k is not too large. 
Combining equations (3.2) and (3.3) we recover the minimal model (1.2). 

The role of the shot noise term P s hot is 1° add power at small length scales, so 
that the ratio P ga \(k) / Pu n (k) is effectively scale-dependent at large k values. This is in 
contrast with the Q model (1.1), which for some values of the nonlinear parameter Q Q \ 
suppresses the power spectrum on the observable scales. Interestingly, taken at face 
value, the P model also predicts a significant /c-dependence for P g3 \(k) / Pn n (k) as 
k — *■ when P s hot once again rises above b 2 Pu n (k) ~ k [29]. It should be noted 
however that this small-A; behaviour is not present for dark matter clustering, since 
momentum conservation demands that the dark matter power spectrum falls off faster 
than P(k) oc k 4 as k — ► [41], a behaviour also observed in numerical simulations [42, 
43]. On the other hand, a non- vanishing shot noise on large scales is in principle not 
forbidden for galaxy clustering, since tracers need not conserve momentum. Whether 
or not this is so remains to be understood. 

In the present work we take the view that since galaxy clustering has not 
been observed on scales where the shot noise behaviour may become problematic 
(k <C 0.01 h Mpc -1 ), equation (1.2) constitutes a sufficient phenomenological model to 
describe the galaxy power spectrum at k > 0.01 h Mpc -1 . 

Additional modifications to the basic P model (1.2) to account for the damping 
of baryon acoustic oscillations and other nonlinear mode coupling effects have been 
discussed in the literature [28,43-45]. These generally lead to nonlinear models of the 
form 

P gal (&) = b 2 A(k)P rm (k) + P shot , (3.4) 

where A(k) is some function of k that depends also on the galaxy type. It is also possible 
to extend the halo model to include redshift space distortion [35,45,46]. We ignore 
these additional corrections in the present analysis in order to keep the number of extra 
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fit parameters to a minimum. However, we emphasise that the effects encapsulated 
by A(k) will become increasingly important as more precise data from future galaxy 
redshift surveys become available. 

4. Test 1: vanilla 

4-1. Set-up 

To compare the performance of the two nonlinear models, we test them against galaxy 
clustering data in three minimal parameter spaces: 

(i) tt m h,tt b h 2 ,h,n s ,ln(l0 10 A), 
(ii) n m h,n b h 2 ,h,n s ,\n(10 w A),P s hot, and 
(hi) Q m h, Q h h 2 , h, n s , \n{10 10 A), Q nl . 

We take the geometry of the universe to be flat, and the initial conditions adiabatic. 
The parameter A = b 2 A s accounts for the normalisation of the galaxy power spectrum, 
and incorporates both the amplitude of the primordial scalar perturbations A s and the 
constant galaxy bias factor b. We do not use any nonlinear correction for parameter 
space (i), i.e., P ga i(k) = b 2 P\\ n {k). For parameter spaces (ii) and (iii), we use the 
P model (1.2) and the Q model (1.1) respectively. 

Note the definition of the matter density parameter Q m h. We choose this 
parameterisation because the turning point of the matter power spectrum is sensitive 
to the comoving Hubble radius at matter-radiation equality, which, for fixed values of 
h, depends on fi m h, not the physical matter density Q m h 2 . For fixed values of n s and h, 
Q m h alone determines the broad shape of the matter power spectrum within the ACDM 
framework. 

We use standard Bayesian inference techniques and the Markov Chain Monte 
Carlo package CosmoMC [47,48] to explore the posterior hypersurfaces as functions 
of the model parameters and the galaxy clustering data sets of section 2. Here, we 
note that within the vanilla framework, the physical baryon density Qt,h 2 , the Hubble 
parameter h, and the scalar spectral index n s can be individually well constrained by 
CMB observations. This information is encapsulated in a set of "WMAP-3 priors" in 
table 1, which we apply when varying these three parameters. This approach differs 
slightly from the more common practice of fixing the parameter values h = 0.72 and 
n s = 1 adopted in, e.g., references [7,26]. Reference [8] further fixes Vt^h 2 = 0.0223. Our 
approach has the advantage that it properly takes into account the uncertainties on these 
parameters and thus avoids two inherent dangers of fixed parameter analyses: biased 
parameter estimates and underestimated errors. Additionally, it permits a consistent 
comparison of our Vt m h constraints not only between different galaxy clustering data 
sets, but also with those obtained from CMB observations alone. 

Broad, top-hat priors are imposed on the remaining parameters (table 1). For 
the parameter Q n \, we choose the range 0-50 for all three data sets. In the absence 
of additional information from, e.g., galaxy formation simulations, the upper limit 
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Table 1. Priors used in our vanilla and vanilla+massive neutrinos analyses. 
Top: WMAP priors based on the three-year data release. These are approximated 
as perfect Gaussians, and we give here their respective mean and standard deviation. 
Bottom: Top- hat priors for the remaining parameters. 



Parameter vanilla vanilla+ncutrinos 

WMAP-3 priors 

n b h 2 0.02229 ±0.00073 0.02158 ± 0.00082 

h 0.732 ±0.032 0.651 ± 0.054 

n s 0.958 ±0.016 0.941 ± 0.022 



Top-hat priors 




1 b m fl 


0.05-0.72 


ln(10 1( U) 


1-5 


-Pshot 


0-min(P o bs) 


Qnl 


0-50 



of this prior is somewhat arbitrary. We will discuss this point in more detail in 
section 6. For the P s hot prior, the upper limit min(P b s ) denotes the minimum clustering 
power measured by a survey. In other words, the linear matter power spectrum 
P\m(k) must not be negative anywhere. For 2dF, SDSS-2 main, and SDSS-4 LRG, 
min(P obs ) = {4000, 1200,9300}, respectively. 

4-2. The internal test: do the individual data sets call for nonlinear correction? 

Figures 1 to 3 show the ID marginal constraints on the parameters {fi m h, P s hot, Qn\} and 
the corresponding minimum \ 2 and numbers of degrees of freedom (d.o.f.) as functions 
of the maximum wavenumber k max included in the analysis. The number of d.o.f. is 
defined here as the number of data points plus the number of priors, minus the number 
of fit parameters; its actual value might be slightly lower, because the power spectrum 
data points are not completely uncorrelated owing to overlapping window functions. 

For fi m h, the green/shaded regions correspond to the ID marginal 68% minimum 
credible intervals (MCI), while the ID modes are indicated by black/dotted lines. On 
the other hand, the marginalised posterior distributions in P s hot and Q n \ are often 
extremely flat, especially at small values of /c max - While it is technically possible to 
construct credible intervals in these cases, quoting such an interval would distract 
from the fact that the parameters are essentially unconstrained. Therefore, instead 
of Bayesian intervals, we give in figures 1 to 3 for P s hot and Q n \ the parameter regions 
satisfying 

-2(lnP-lnP max )<l, (4.1) 

where V is the ID marginal posterior and P max denotes its value at the ID mode. We 
loosely label this the "la" interval. For large A; max values, especially k max ~ 0.2 h Mpc -1 , 
where the posterior distributions are approximately Gaussian, this la interval is identical 
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Figure 1. ID marginal 68% MCIs for fl m h and la intervals (see main text for 
definitions) for the nonlinear parameters P s hot and Q n \ as functions of A: max using the 
SDSS-4 LRG power spectrum (green/shaded regions). Left: No nonlinear correction. 
Centre: Correction with the P model. Right: Correction with the Q model. In each 
case the black dotted line indicates the ID mode. The top plot in each column 
shows the minimum \ 2 (red/solid) and the number of degrees of freedom in the fit 
(blue/dash). 



to the 68% MCIs. See reference [49] for more detailed discussions of the various 
statistical quantities. 

Changes in the Vt m h estimates The effects of nonlinear correction are most evident 
when we include data beyond k max ~ 0.1 h Mpc -1 . Here, nonlinear correction generally 
shifts the Q m h estimates to lower values relative to the case with no correction, 
irrespective of the nonlinear model used. For SDSS-4 LRG, this shift is very dramatic: at 
kmax ~ 0.2 h Mpc -1 the 68% MCI moves down by an amount comparable to six or seven 
times its half-width (figure 1). In contrast, the shifts induced for 2dF and SDSS-2 main 
by either nonlinear model are mild: at fc max ~ 0.2 h Mpc -1 a small overlap between the 
68% MCIs before and after correction can still be seen in figures 2 and 3. 



The x 2 t es t The need for nonlinear correction in the case of SDSS-4 LRG is 
corroborated by a comparison of the minimum \ 2 values. Between correction and no 
correction, figure 1 shows that at /c max ~ 0.2 h Mpc -1 , Xmin changes from ~ 60 for 
18 d.o.f. to ~ 20 for 17 d.o.f., again irrespective of the nonlinear model used. On the 
other hand, similar comparisons for 2dF and for SDSS-2 main in figures 2 and 3 do not 
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Figure 2. Same as figure 1, but for 2dF. 
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Figure 3. Same as figure 1, but for SDSS-2 main. 
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indicate any urgent need for an extra correction parameter: for 2dF, A%^ in ~ 4; for 
SDSS-2 main, A%^ in is virtually negligible. 



Internal consistency We consider a data set internally consistent if the Q m h credible 
intervals agree for all choices of /c max . For SDSS-4 LRG, figure 1 clearly demonstrates 
that nonlinear correction is necessary in order to achieve some semblance of internal 
consistency. A small discrepancy remains after correction: at /c max ~ 0.03 and 
0.2 h Mpc^ 1 the 68% MCIs do not quite touch each other. However, this discrepancy 
is not statistically significant. We bring up this point here nonetheless as a cautionary 
note against over-interpretion of credible intervals. 

Nonlinear correction also helps to improve consistency in the Q m h estimates from 
2dF in figure 2. As for SDSS- 2 main, the large uncertainties in Q m h in figure 3 means 
that internal consistency is not an issue, with or without nonlinear correction. 

How much correction? Focussing on the cases with nonlinear correction, we see in 
figures 1 to 3 that while SDSS-4 LRG clearly prefers a nonzero correction parameter — 
be it -Pghot or Qui — at A; max > 0.1 h Mpc -1 , the 2dF and the SDSS-2 main data sets do 
not exhibit the same strong preference. For SDSS-2 main and most choices of /c max for 
2dF, the la regions include P s hot = and Q n \ = (although, as mentioned before, the 
case of Q n \ = is not equivalent to no nonlinear correction). At A: max < 0.1 h Mpc -1 , 
the parameters P s hot and Q Q \ cannot be constrained by data. 

In terms of the P model, roughly half of the small scale (k > 0.1 h Mpc -1 ) power 
in SDSS-4 LRG can be attributed to the shot noise term. For 2dF and SDSS-2 main, 
the shot noise contribution is about a quarter according to the ID mode. 

Which nonlinear model? Perhaps the most striking feature about the two nonlinear 
models considered here is that their corrective effects appear to be identical. This is 
particularly apparent in figures 1 and 2 at A; max > 0.1 h Mpc -1 , where the preferred 
values of the correction parameters P s hot and Q n \ exhibit virtually the same dependence 
on A; max . At fc max ~ 0.2 h Mpc -1 the VL m h estimates also show little if any dependence 
on the priors imposed on the nonlinear correction parameters. Thus as far as vanilla 
cosmology is concerned, there is no preference for either nonlinear model from a 
phenomenological standpoint, although the transparency of the P model still makes 
it the more attractive one of the two. 

4-3. The external test: are all data sets consistent with each other? 

Figure 4 shows how the constraints on Q m h obtained from different galaxy clustering 
data sets and with different nonlinear correction methods compare with each other. 
For good comparison we indicate in the figure also the corresponding estimate from 
WMAP-3. Similar information is available in table 2, in which we give the ID marginal 
68% and 95% MCIs at k nmx ~ 0.2 h Mpc" 1 (fc max ~ 0.18 h Mpc" 1 for 2dF). 
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Figure 4. ID marginal 68% MCIs for £l m h as functions of fc ma x using galaxy clustering 
data from SDSS-4 LRG (red/shaded region), 2dF (blue/dashed line) and SDSS-2 main 
(green/solid line). Top: No nonlinear correction is used. Middle: Correction with the 
P model. Bottom: Correction with the Q model. For comparison we show also the 
corresponding 68% MCI from WMAP-3 (grey/dotted line). 



For both SDSS-4 LRG and SDSS-2 main, nonlinear correction clearly leads to 
better agreement with 2dF and with WMAP-3 at /c max > 0.1 /i Mpc -1 . In the case 
of SDSS-2 main, although the correction is not sufficient to cause the 68% region to 
overlap with that from 2dF (using the same correction method), the two 95% regions 
are certainly consistent. Importantly, the level of disagreement between SDSS-2 main 
and 2dF after correction is no worse than the small internal inconsistency between 
the low and the high fc max constraints on Q m h from SDSS-4 LRG already discussed in 
section 4.2. Thus if we should accept that the nonlinear models (1.1) and (1.2) offer 
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Table 2. ID marginal 68% (95%) MCIs for Sl m h and the nonlinear parameters P s hot 
and Q n i from the three galaxy clustering data sets under consideration. We show 
results for fc max ~ 0.2 h Mpc -1 (fc ma x ~ 0.18 h Mpc -1 for 2dF), using (i) no nonlinear 
correction, (ii) correction with the P model, and (iii) correction with the Q model. 
The WMAP-3 preferred regions are also quoted here for comparison. 



Data set 


No correction 


P model 


Q model 


1 lyjltl 

SDSS-4 LRG 


0.234-0.252 (0.227-0.262) 


0.163-0.184 (0.153-0.195) 


0.169-0.191 (0.158-0.202) 


2dF 


0.167-0.204(0.152-0.228) 


0.134-0.173 (0.122-0.196) 


0.134-0.178 (0.118-0.207) 


SDSS-2 main 


0.233-0.277 (0.213-0.300) 


0.190-0.250 (0.167-0.279) 


0.193-0.276 (0.167-0.316) 


WMAP-3 




0.159-0.195 (0.142-0.211) 




Phot or Q n l 








SDSS-4 LRG 


- 


3980-5170 (3380-5740) 


22.7-25.5 (18.6-37.9) 


2dF 


- 


452-1310 (< 1640) 


7.0-19.2 (2.5-27.3) 


SDSS-2 main 


- 


< 506 (< 831) 


< 10.2 (< 19.5) 



sufficient correction for SDSS-4 LRG, logically we must also consider them adequate for 
SDSS-2 main. 

Lastly, we observe in figure 4 that the 2dF preferred values of Q m h tend in any 
case to be on the low side of SDSS-2 main, even at values of /c max well below those at 
which nonlinearity nominally sets in. This suggests that the residual inconsistency after 
nonlinear correction can rather be put down to a statistical aberration at small k values, 
than is indicative of a failure of either nonlinear model. Indeed, Sanchez and Cole [26] 
compared directly the power spectra of red galaxies (believed to be the main source of 
scale-dependent biasing) from 2dF and the main galaxy sample of SDSS data release 5, 
and found a similar discrepancy in the raw power spectrum data (see figure 7 of their 
paper). Our SDSS-2 main data set is but a subsample of SDSS data release 5; that it 
contains the same small fluctuation should be of no surprise. 

5. Test 2: vanilla+massive neutrinos 

It is well known that a subdominant component of massive neutrino HDM in the 
matter content slows down the growth of density perturbations on small length scales. 
In terms of the matter power spectrum, we expect a suppression of power of order 
APi in /Pii n ~ 8 Q u /Q m at k ^> A; F s, where fcps is the neutrinos' free-streaming 
wavenumber, and Q u the neutrino energy density. Since this suppression effectively 
changes the shape of the matter power spectrum at large wavenumbers, some amount 
of degeneracy could conceivably exist between the neutrino energy density and the 
nonlinear correction parameters. In this section we investigate the possible existence of 
such a degeneracy, and if so, its effects on massive neutrino cosmology. 



Nonlinear corrections to the cosmological matter power spectrum 13 

Table 3. ID marginal 68% (95%) MCIs for Q m h, /„, and the nonlinear parameters 
Pshot and Q n i from SDSS-4 LRG, and the associated minimum \ 2 values. These 
are obtained using (i) no nonlinear correction, (ii) correction with the P model, and 
(iii) correction with the Q model. We quote here also the WMAP-3 preferred regions 
for comparison. 



Model Cl m h f v P s hot or Q n \ xl 



nun 



No correction 0.246-0.287 (0.232-0.314) 0.01-0.047 (< 0.065) 63.2 

P model 0.168-0.216 (0.154-0.256) < 0.057 (< 0.112) 4070-5270 (3430-5810) 19.7 

Q model 0.171-0.226 (0.158-0.274) < 0.061 (< 0.135) 23.6-36.9 (17.5-47.5) 20.6 



WMAP-3 0.178-0.245 (0.151-0.270) < 0.080 (< 0.127) 



5.1. Set-up 

We consider three parameter spaces: 

(i) f u ,n m h,n b h 2 ,h,n s ,\n(10 10 A), 
(ii) f u ,D, m h,D, b h 2 ,h,n s ,ha.(10 w A), P shot , and 
(iii) /„, Q m h, Q b h 2 , h, n s , m(10 1( U), Q nl . 

Here, the neutrino fraction f v = Q u /Q m is defined as the ratio of the neutrino energy 
density Q u to the total matter density Q m . The former is given by the well known 
expression 

where J2 m v denotes the sum of the neutrino masses. We assume 3.04 degenerate 
neutrino species, which should be a good approximation since neither CMB nor galaxy 
clustering measurements are at present sufficiently sensitive to X)m„ ;$ 0.3 eV, where 
one might expect some small effects due to the neutrino mass hierarchy. 

We fit these three cases to SDSS-4 LRG up to k max ~ 0.2 h Mpc -1 , using no 
nonlinear correction in case (i), and correction with the P and the Q model respectively 
for cases (ii) and (iii). As in the previous section, we impose WMAP-3 priors on the 
parameters h, n s , and Q b h 2 tabulated in table 1. Note that these priors differ from those 
used previously, since CMB constraints are model-dependent. We adopt a top-hat prior 
on /„, 0-0.5, while for Q n \ we use 0-100, in anticipation that more nonlinear correction 
may be required to offset the higher Q m h values usually inferred in cosmologies with 
massive neutrinos. 

5.2. Results and discussions 

Table 3 shows the ID 68% and 95% MCIs for Q m h, /„, and the nonlinear parameters 
P s hot and Qni for the three cases considered. We also give the minimum x 2 values as a 
measure of the goodness-of-fit. 
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Figure 5. 2D marginal 68% and 95% MCIs for {/„,P sho t} and {/„, Q n i} from 
SDSS-4 LRG, for cases (ii) and (iii) of section 5.1. 



As in the case for vanilla cosmology, fitting the SDSS-4 LRG data without nonlinear 
correction leads to a very poor goodness-of-fit; the Xmin varue is 63.2, for approximately 
20 + 3 — 6 = 17 degrees of freedom (20 data points of SDSS-4 LRG, 3 for the priors 
on n s , h and Qf,h 2 , and —6 for 6 free parameters). Adding nonlinear correction reduces 
Xmin by more than 40 units at the expense of only one extra parameter. Again, both 
nonlinear models offer very similar corrections to the power spectrum in terms of the 
inferred Q m h and f v values, although the P model appears to provide a slightly better 
fit to the data judging by its slightly smaller Xmin- Interestingly, despite the dramatic 
decrease in the minimum x 2 between correction and no correction, the resulting shifts in 
the Q m h and/y estimates are deceptively small so that the 95% MCIs remain compatible 
before and after correction. 

In the absence of nonlinear correction, the 95% upper limit on f u is about a factor 
of two too tight compared with the corrected case. This situation is reminiscent of 
the overly constraining bounds on J2 m v derived in some recent combined analyses of 
WMAP-3 and the flux power spectrum of the Lyman-a forest [50]. Too much power 
at large wavenumbers — either because of uncorrected nonlinearities in the galaxy power 
spectrum or an unusually large normalisation in the case of the Lyman-a forest — leads 
to the appearance of an overly flat matter power spectrum, which in turn prefers a 
smaller neutrino fraction and hence a smaller neutrino mass. The difference between 
the two corrected f v bounds is about 20%. The corresponding 95% limits on Yli m v-, 
< 1.76 eV for the P model and < 1.83 eV for the Q model, differ by even less. 
Thus cosmological neutrino mass determination is at present unaffected by our choice 
of nonlinear correction model. 

Finally, we see in figure 5 that no strong degeneracy exists between the neutrino 
fraction f v and the nonlinear correction parameters P s hot and Q n \- For small values 
of f v , the changes induced in the linear matter power spectrum by massive neutrino 
dark matter can be approximately mimicked by a redefinition of the apparent Vt m h 
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parameter [7], 

(,' 'm^Japparcnt (,' 'm"ytrue *--^Jvi V / 

assuming h ~ 0.7. Here, (^m^)apparcnt is the parameter that is actually constrained by 
power spectrum data, while (fl m /j)t rue denotes the true value of Q m h. This expression 
also encapsulates the well known (approximate) degeneracy between Q m h and /„. Larger 
f v values induce more complicated changes in the power spectrum, and it is reassuring to 
see that these changes are not degenerate with nonlinear correction using either model. 

6. The pathology of the Q model: vanilla+thermal axions 

While both nonlinear models work very well for vanilla and for vanilla+massive neutrino 
cosmologies, and we may be tempted to conclude that they are for all purposes 
phenomenologically identical, we provide in this section a counter-example in which 
the incorrect functional form of the Q model can lead to some misleading results. 

The case in point is a class of models containing a possible subdominant HDM 
component due to relic thermal axions with mass m a . These models differ from those 
with massive neutrino HDM in that the temperature T a = [\0.7h/ g* D {m a )] 1 ^ T v and 
hence the number density n a = [C(3)/7r 2 ] T 3 of the axions are functions of m a , since 
m a determines when the particle species should decouple from the primordial plasma. 
Here, the function g* D (m a ) denotes the effective number of thermal degrees of freedom 
at the time of decoupling, and must be calculated by carefully tracking the freeze-out 
process (e.g., [16]). 

Thus, although qualitatively thermal axion HDM exhibits free-streaming features 
very similar to those of massive neutrinos, quantitatively the suppression of small scale 
power in the matter power spectrum has a nonlinear dependence on the axion mass. 
We show in this section how this nontrivial dependence can cause some problems for 
the Q model. 

6.1. Set-up 

We consider three parameter spaces: 

(i) m a , Q m h, fibh 2 , h, n s , ln(10 10 A,), b, Q n \, with a top-hat prior 0-100 on Q nh 
(ii) m a , Q m h, fibh 2 , h, n s , ln(10 10 A s ), b, Q n \, with a top-hat prior 0-200 on Q nh and 
(hi) m a , n m h, Q b h 2 , h, n s , ln(10 10 A s ), b, P shot . 

Cases (i) and (ii) both use the nonlinear model (1.1), and differ only in the priors 
imposed on the correction parameter Q n \. Case (iii) uses the nonlinear model (1.2), 
with the usual prior on P s hot (see table 1). We fit each case to the combined data set 
WMAP-3+SDSS-4 LRG, using data up to k max ~ 0.2 h Mpc" 1 for the latter. 
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Figure 6. 2D marginal 68% and 95% MCIs for {il m h, m a }, {£l m h, Qni}, and {m a , Qni} 
from WMAP-3+SDSS-4 LRG. In the diagonal are the ID marginal posteriors for the 
three named parameters. Black/solid lines correspond to case (i) of section 6.1 with 
a top-hat prior 0-100 on Q n \, while the blue/shaded regions and blue/dashed lines 
denote case (ii) with Q n i prior 0-200. 



6.2. Results and discussions 

Figure 6 shows the 2D marginal 68% and 95% MCIs for {Q m h,m a }, {Q m h,Q n \}, and 
{m a ,Q n \}, as well as the ID marginal posteriors for the the same three parameters for 
cases (i) and (ii). 

Consider first the ID marginal posteriors for Q n \. We see in figure 6 that the 
posteriors in both cases (i) and (ii) remain finite all the way up to the upper limit of 
the prior imposed on Q n \. This is also reflected in the relevant 2D contours, which are 
abruptly cut off at Q n \ = 100 and Q n \ = 200 respectively Data alone does not constrain 
Qni in this class of cosmological models. 

While this does not affect the Q m h estimates, the inference of the axion mass m a 
depends crucially on how well we can constrain Q n \ because of a persisting degeneracy 
between the two parameters. Indeed, if we do not cut off Q n \ by hand at some sufficiently 
small value, a second peak begins to appear in the posterior at {m a ~ 6 eV, Q n \ ~ 200}, 
besides the one at {m a ~ eV, Qni ~ 30}. From figure 7 we see that the power spectra 
for the two peaks after nonlinear correction are almost perfectly degenerate up to (and 
beyond) k ~ 0.2 h Mpc^ 1 , even though the corresponding linear matter power spectra 



Nonlinear corrections to the cosmological matter power spectrum 

,5 



17 



^ 



D. 10 



-4£ 








10 



0, 10 



5 



4 



3 



2 



10" r 



10 



1 



' ' ' ' ' ""' 




r ^-Corrected 


1 


! m a - 0.6 eV 


// 


I ■ m a = 5.9 eV 


i 





~ 



0.001 0.010 0.100 1.000 
k(h/Mpc) 

Figure 7. Matter power spectra for the two peak values of m - The black/solid lines 
are for {m a = 0.6 eV, Q n \ — 31}, and the red/dashed lines for {m a = 5.9 eV, Q n \ = 
170}. The upper panel shows the linear matter power spectra, while the lower 
panel includes nonlinear correction with the Q model. The vertical line indicates 
the maximum value of k used in the analyses. 



differ markedly already at k ~ 0.04 h Mpc -1 . 

We may be inclined to regard a 6 eV axion and, by implication, a large Q n \ 
value as unphysical because the former runs in conflict with constraints on m a derived 
from stellar energy loss arguments and from telescope searches for axion radiative 
decays [51]. Furthermore, the m a , Q n i-degeneracy exists only in a limited k range: 
as shown in figure 7, the corrected power spectra for the two peaks begin to deviate 
at k > 0.3 h Mpc -1 . This suggests that the exact location of the high m a peak in 
{m a , Q n i}-space may in fact depend on our choice of fc max - Naturally we could have 
avoided this second peak simply by demanding consistency with other astrophysical 
constraints on m a . However, in the absence of, e.g., galaxy formation simulations for 
this very class of cosmological models, we have a priori no reason to reject Q n \ ~ 200 or 
any other higher or lower value besides our own prejudices. Moreover, even if we manage 
to avoid the second peak with a finely tuned prior on Q nh the m a , Q n i-degeneracy means 
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Figure 8. 2D marginal 68% and 95% MCIs for {Q m h,m a }, {n m h,P s hot}, and 
{m a ,P shot } from WMAP-3+SDSS-4 LRG for case (iii) of section 6.1. ID marginal 
posteriors for the same three parameters are shown in the diagonal. 



that the constraint thus obtained on m a will still depend sensitively on exactly what we 
choose to be the upper limit of that prior. 

This exercise also highlights the danger of analytic marginalisation [47], a technique 
sometimes used on nuisance parameters in CosmoMC to shorten the computation time. 
Here, analytic integration of the posterior in the direction of a nuisance parameter is 
made possible by taking the lower and upper limits of the parameter's top-hat prior 
to — oo and oo respectively. Using this technique on the marginalisation of Q n \, we 
find a unimodal ID posterior for m a whose 95% MCI of 4.95 < m a /eV < 7.17 favours 
unambiguously the high m a region. Thus, analytic marginalisation can be very useful if 
the likelihood function itself sufficiently constrains the nuisance parameters. Otherwise, 
as we have seen here, the end results are potentially misleading. 

Finally, we note that the P model does not suffer these problems. Figure 8 shows 
the 2D marginal 68% and 95% MCIs for {Q m h, m a }, {Q m h, -Pshot}, and {m a , P s hot}, and 
the corresponding ID marginal posteriors for case (iii). Here, although the correction 
parameter P s hot is slightly degenerate with the axion mass m a , it is independently well 
constrained by data. Importantly, even if data fails to constrain P s hot 5 we have a simple 
and well defined way to choose our prior on P s hot- Thus we conclude that the P model is 
at present superior to the Q model for nonlinear correction in non-vanilla cosmologies. 
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7. Conclusions 

In this paper we have explored two nonlinear correction and scale-dependent bias 
models — the Q model of reference [7] and the halo model-inspired P model — in some 
detail. We have confronted them with a range of galaxy clustering data and cosmologies 
to determine the strengths and weaknesses of the models. 

In the context of standard ACDM cosmology, we find that both models perform 
equally well on present galaxy power spectrum data, in the sense that their corrective 
effects on the matter power spectrum are essentially identical. The case for nonlinear 
correction is very strong for the SDSS-4 LRG power spectrum. An indicative figure 
of merit is the minimum x 2 : fitting data up to A: max ~ 0.2 h Mpc -1 , we find that 
Xmin changes from ~ 60 for 18 degrees of freedom without correction to ~ 20 for 17 
degrees of freedom with correction, irrespective of the nonlinear model used. For 2dF 
and SDSS-2 main, however, the need for correction is marginal and subsists primarily 
because their respective Q m h estimates show better agreement with than without 
nonlinear correction. The preferred values of Q m h from SDSS-4 LRG, SDSS-2 main, and 
2dF after correction, as well as from WMAP-3 can all be reconciled at 95% confidence, 
contrary to the case without correction. 

Similar results are obtained for ACDM cosmologies extended with a subdominant 
component of massive neutrino hot dark matter. Data again show no strong preference 
for either nonlinear correction model, nor do we find any detrimental degeneracy between 
the neutrino fraction /„ and either nonlinear correction parameter. Cosmological 
neutrino mass determination is at the time being unaffected by our choice of nonlinear 
correction model. 

However, if the subdominant free-streaming dark matter is in the form of relic 
thermal axions, a nontrivial degeneracy between the axion mass m a and the correction 
parameter Q n \ renders the Q model highly pathological, so that our inference of m a 
depends sensitively on the prior we impose on Q nl . In contrast, the P model, whose 
functional form is based on well motivated physics, does not suffer from this problem, 
and is arguably superior to the Q model. Note that we have used relic thermal axions 
here as an example. But our results may also be relevant for other light thermal relics 
whose temperature and hence abundance depend on the mass of the particle. 

More precise data from future galaxy redshift surveys will eventually render 
these simplistic models inadequate to describe nonlinear evolution and scale-dependent 
biasing. For example, the damping of baryon acoustic oscillations due to nonlinear 
mode-coupling will need to be factored into the game at some stage [28,43,44]. The 
advent of wide- and deep- field weak gravitational lensing surveys in the next decade as an 
alternative probe of the large scale structure distribution will circumvent some of these 
nonlinearity issues. But galaxy redshift surveys will remain an important tool for the 
observation of baryon acoustic oscillations, and nonlinear evolution/scale-dependent bias 
modelling will continue to constitute an important aspect of the cosmological analysis 
machinery. 
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